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^ 1 ABSTRACT 
> j 

A classical problem of free-space Green's function C7oa representations of the Helmholtz 
equation is studied in various quasi-periodic cases, i.e., when an underlying periodicity is 
imposed in less dimensions than is the dimension of an embedding space. Exponentially 
convergent series for the free-space quasi-periodic C? a and for the expansion coefficients 
Dl of C7 a in the basis of regular (cylindrical in two dimensions and spherical in three 
. dimension (3D)) waves, or lattice sums, are reviewed and new results for the case of a 

^! one-dimensional (ID) periodicity in 3D are derived. From a mathematical point of view, 

a derivation of exponentially convergent representations for Schlomilch series of cylindri- 
cal and spherical Hankel functions of any integer order is accomplished. Exponentially 
convergent series for Goa and lattice sums Dl hold for any value of the Bloch momentum 

>• ■ and allow C7 a to be efficiently evaluated also in the periodicity plane. The quasi-periodic 

Green's functions of the Laplace equation are obtained from the corresponding represen- 
tations of Goa of the Helmholtz equation by taking the limit of the wave vector magnitude 
going to zero. The derivation of relevant results in the case of a ID periodicity in 3D high- 
lights the common part which is universally applicable to any of remaining quasi-periodic 
cases. The results obtained can be useful for numerical solution of boundary integral 
equations for potential flows in fluid mechanics, remote sensing of periodic surfaces, peri- 
odic gratings, and infinite arrays of resonators coupled to a waveguide, in many contexts 
of simulating systems of charged particles, in molecular dynamics, for the description of 
quasi-periodic arrays of point interactions in quantum mechanics, and in various ab-initio 
first-principle multiple-scattering theories for the analysis of diffraction of classical and 
quantum waves. 

PACS numbers: 02.70.Pt; 03.65.Db; 03.65.Nk; 34.20.-b; 42.25.Fx; 46.40.Cd; 47.15.km; 
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1 Introduction 



Let A be a ^A-dimensional simple (Bravais) periodic lattice embedded in a space of di- 
mension d > d\ (the condition of a simple lattice can easily be relaxed to an arbitrary 
periodic lattice by following recipes of Refs. [H El El IH ISj)- Let Hq stand for the cylin- 
drical (Hq 1 ^) and spherical (h^) Hankel functions in d = 2 and d = 3, respectively, be 
corresponding angular momentum harmonics (cylindrical in d — 2 and spherical in d = 3; 
see also Appendix \B§ for properties of 3^l) , and r and r' be spatial points. The article is 
concerned with an efficient calculation of the series 

Y,nt(a\r-r' + r n \)yi(r n )^, (1) 

r„6A 

where the origin of coordinates is in the lattice, k is called the Bloch momentum, a = 27r/ A 
denotes a wave vector magnitude (a is not necessarily equal to |k|), A is a wavelength, 
and L is, in general, a multi-index of angular momentum numbers [e.g., L = (Im) in three 
dimensions with I > and —l<m<l\. 

In mathematical literature such a series are known as Schlomilch series [HI El IE] • As 
first noted by Emersleben [0] , in three dimensions (3D) in the special case of r — r' = and 
/ = a = the series reduce to Epstein zeta functions jlUt fTT| fT^| . A physical motivation to 
investigate such a series derives from a fact that, for / = 0, the series (0) are prerequisite 
to determine a corresponding free-space (quasi-) periodic Green's function G a of a scalar 
Helmholtz equation (see below). For I ^ and r — r' = 0, with singular term being 
excluded, the series then formally determine the lattice sums Dl, which are defined 
as the expansion coefficients Goa i n the basis of cylindrical and spherical Bessel functions 
in 2D and 3D, respectively [see Eqs. (JSJ), (fTTJl below]. 

However, analytic closed expressions for such series are only known for I = in two 
particular cases in 3D: (i) in the case of a one-dimensional (ID) periodicity [THJ [T^l [TH] and, 
when additionally a = (the Laplace limit), (ii) for a two-dimensional (2D) periodicity 
|14j . Otherwise the summation in (Q) has to be performed numerically. However, since 

/n\(<W)/J 

Ht{z) ~c h z- {d - 1)/2 f-J exp{z[^- {d- 1)tt/4 - Ztt/2]} (2) 

as z — > oo and — ir < argz < 2tt, where Ch = 1 unless d = 2, in which case = y/2/ir (see 
Eqs. (9.2.3) and (10.1.1) of Ref. [E]), the series (JTJ) is not absolutely convergent. Even 
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if one assumes that a has an infinitesimally small positive imaginary part, and thereby 
establishing absolute convergence, the convergence of the series in Eq. is notoriously 
slow, thereby rendering it useless for practical applications. 

The study of efficient techniques for the calculation of G a and lattice sums D L has 
a long history [U El S El 113 dl EH EDI ED 1221 123 El and the topic has roots and 
branches in widely differing areas of chemistry, physics, and mathematics. Despite that 
it still continues to be a perennial research subject 0E3EHlEllEaEHlllIlEIlE21EIIl 
EH EH EH EZI- However, quite often this happens merely because researches in widely 
differing areas of chemistry, physics, and mathematics are not aware of the techniques 
and results developed in connection with the so-called Korringa-Kohn-Rostoker (KKR) 
PU 121 123 121] and layer Korringa-Kohn-Rostoker (LKKR) theories either for quantum 
(electron) waves within low-energy electron diffraction (LEED) theory [1J EDI El E2] or 
for various classical (acoustic, elastic, electromagnetic, water) waves [SlEniEZl- Therein 
exponentially convergent series for Goa and lattice sums Dl have been derived for a number 
of cases. In a periodic case (g?a = d), efficient computational schemes for Goa and lattice 
sums Dl have been provided more than thirty years ago [21 ESI- The quasi-periodic case 
g?a = 2 and d = 3 has been investigated in detail in a series of articles by Kambe almost 
forty years ago [UEII1E]- The case d^ = 1 and d = 2 has been dealt with only relatively 
recently in Refs. j3 ESI- Surprisingly enough, the case g?a = 1 and d = 3 has not 
been studied in full detail yet, although it may provide an efficient description of "wave 
integrated circuits" . 

Therefore, in the following we shall focus on the so-called quasi-periodic, or layer, 
case, when the underlying lattice A is of lower dimensionality than the embedding space 
(g?a < d), and, in particular, on the d\ — 1, d = 3 case. There are many physical 
problem which would profit from an efficient computational scheme for Goa and lattice 
sums Dl in the case of ID periodicity in 3D. For instance, a Green's function representing a 
point source and satisfying the respective von Neumann and Dirichlet boundary condition 
on a flow channel walls in fluid mechanics for the flow between parallel planes can be 
written as a sum and difference of Goa corresponding to ID periodicity in 3D taken at 
two different spatial points (SHI- (For the flow in a rectangular channel, the relevant 
Gqa would then correspond to 2D periodicity in 3D [SHI-) Another problem involves 
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infinite arrays of resonators coupled to a waveguide which by itself may have a plethora of 
photonics applications [3B] • Moreover, recent progress in nanotechnology made it possible 
to fabricate ID chains of metal nanoparticles [331101 and dielectric microparticles |41j . 
Additionally, understanding of linear periodic arrays of lossless spheres is germane for 
the qualitative description of a finite- lengths periodic arrays of small antennas [32|. In a 
linear chain of spherical metal nanoparticles light can be transmitted by electrodynamic 
interparticle coupling resulting in a subwavelength-sized light guide [3*"*"1 f4T"l 14*31 1^1- 
So far, infinite ID linear chains of particles have only been investigated within the one 
particle theory framework of Schrodinger equation for the description of polymers [T!*"j, in 
the electrostatic limit within the framework of the Laplace equation [4*4*| . or in the dipole 
approximation [4*2*1 ["""""J [!"""] . Moreover, the energy operator in the one-particle theory of 
periodic point (zero-range) interactions is constructed in terms of an auxiliary operator 
which corresponds essentially to the operator of multiplication by D 00 (7 function of 
Karpeshina [T3*l I47j ) . 

In the following, exponentially convergent series for the free-space quasi-periodic Goa 
and for the lattice sums Dl are reviewed and new results for the case of a ID periodicity in 
3D are derived [see Eqs. (|83|). (jlU2j) . and (jll8|) below]. The derivation of relevant results 
for the case of a ID periodicity in 3D is performed in such a way that the common part 
which is universally applicable to any of the remaining quasi-periodic cases is highlighted. 
Thereby a link with earlier results by Kambe ""fl [201 f° r a 2D periodicity in 3D can easily 
be established and the proof of results for a ID periodicity in 2D announced in Ref. [213 
can easily be carried out. 

1.1 The outline of the article 

The article is organized as follows. Sec. [21 introduces notation, provides some necessary 
definitions and gives an overview of some of the problems requiring the knowledge Goa 
and lattice sums. In Sec. El Goa is expressed as an exponentially convergent sum over a 
dual lattice A*. Such a dual representation of Goa is then a starting point in the derivation 
of a corresponding exponentially convergent Ewald representation of Goa i n Sec. [3] As 
in the bulk derivation of the Ewald representation invokes a suitable integral 

representation of Hankel functions and a Jacobi identity (see Appendix \D\ below) . In 
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addition, following Kambe [21], an analytic continuation procedure is employed, which is 
analogous to finding the values of the Riemann ^-function outside the domain of absolute 
convergence of a defining series (see Ref. [IB], p. 273). The Ewald representation of G a, 
which is a hybrid sum over both A and dual lattice A* , converges uniformly and absolutely 
over bounded sets of R. Unlike the dual representation, the Ewald representation can also 
be efficiently evaluated in the periodicity plane (provided that it remains R ^ A). 

Exponentially convergent series for the lattice sums in the case of a ID periodicity 
in 3D are derived in Sec. El Eqs. (jHSJ), (jl(J2j) . ()118|) of the section are the main new results 
of this paper. 

In Sec. El the quasi-periodic Green's function of the Laplace equation are obtained 
from that of the Helmholtz equation by taking the limit o —>■ 0. We then end up with 
discussion (Sec. |7J) and summary and conclusions (Sec. |SJ). 

To make this paper as readable as possible, several technical arguments have been 
relegated to a number of appendices. Appendix 1X1 summarizes relevant integral represen- 
tations of Til and Appendix [B] lists relevant properties of harmonics 3^l- Some useful 
properties of free-space scattering Green's function are collected in Appendix whereas 
Appendix iDl shows several forms of Jacobi identities. In Appendix [E] the dual representa- 
tion of the quasi-periodic Green's functions is derived by following original path of McRae 
[4U] for a 2D periodicity in 3D, i.e., by applying an Ewald integral representation and the 
generalized Jacobi identity. Some of the general properties of free-space quasi-periodic 
Greens functions and the lattice sums are outlined in Appendix [Fj Alternative definitions 
of lattice sums and structure constants are then summarized in Appendix [UJ 

2 Notation and definitions: Goa and lattice sums 

2.1 G A 

A corresponding free-space (quasi-) periodic Green's function G a of a scalar Helmholtz 
equation is defined by an image-like series 

G 0A (a, k, R) = R + r„)e- lk - r " = £ G+{a, R - r„)e lk ' r ", (3) 

r„eA r„eA 
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where the origin of coordinates is in the lattice. Here Gq is the free-space scattering, or 
retarded, Green's function of a scalar Helmholtz equation in d dimensions, 

[A + a 2 ]G (a,r,r')=5(r-r'), (4) 

A being a corresponding Laplace operator, which is represented for large R = |R| = |r— r'| 
by outgoing waves (i.e. satisfies Sommerfeld radiation condition). Since Gq is only a 
function of R = r — r', the functional dependence of Gq has been written as Gq(u, R) in 
Eq. ©. 

Gq is proportional to an appropriate Hankel function of zero order ^E] 

G+(a, R) = lim — i-r f „ e t R dk=-£ — ^ a^HUaR) (5) 

(see also Appendix EJ) • Here k = |k| and A is the surface of a unit sphere in <i-dimensions 
[see Eq. (I147|) below]. (With TCq(z) = e * z as a one- dimensional (ID) analog of the Hankel 
function [SD] Eq. (JHJ) becomes also valid for d = 1, in which case A = 2 (unit "sphere" in 
ID consists of two points).) Consequently, upon substituting (jHJ) into (jHJ one arrives up 
to a proportionality factor to the special case of series (JTJ for / = 0. 

A scalar Helmholtz equation is employed for the description of various waves arising 
in acoustics, mechanics, fluid dynamics, electromagnetism, and quantum mechanics [51J. 
An important class of problems which requires an efficient calculation of Gq\ arises in 
connection with remote sensing of periodic surfaces [23], numerical solution of boundary 
integral equations for potential flows in fluid mechanics [2E1 EH EDI IHJ] , periodic gratings, 
and infinite arrays of resonators coupled to a waveguide jHH], i n descriptions of dipolar 
fields in simulated liquid- vapour interfaces in many contexts of simulating systems 
of charged particles, such as crystal binding and lattice vibrations [3 HI], Madelung 
constant [TT) IT2"]. in various problems in molecular dynamics and Monte Carlo simulations 
of particles interacting by long-range Coulomb forces jH2 EH] , wherein periodic boundary 
conditions are usually imposed in order to avoid the boundary effects. 

2.2 Lattice sums 

Within a primitive cell of A, the respective Green's functions Gq and Goa only differ up 
to boundary conditions and their respective singular parts are identical. For the scalar 



5 



Helmholtz equation, the singular, or principal- value, part of Gq is 

GS (,,R) = Re G +(,,R) = { _ co ^/ ( t R) , Id < 6 > 

where Nq is the cylindrical Neumann function ^H]- Therefore, the difference 

D A (a, k, R) = Gq\(<j, k, R) - G p (a, R) (7) 

is regular for R — ► and can be expanded in terms of regular [cylindrical in 2D, spherical 
in 3D] waves H H3 EH ESI E3 E3 

Av(<7,k,R) = ^2D L (a,k)J l (aR)y L (R). (8) 

L 

Here the symbol J\ stands for cylindrical and spherical Bessel functions in 2D and 3D, 
respectively. The expansion coefficients Dl(ct, k) introduced by Eq. (JHJ) are then the 
sought lattice sums. The choice of what to subtract off in Eq. (JJJ) is somewhat arbitrary 
and other choices will lead to slight amended expressions for the lattice sums Di (see 
also Appendix |GJ) • The present choice goes back to Kohn and Rostoker [33] and has been 
adopted by Ham and Segall 0, Kambe [31 H3, Pendry [22], and others [2H1I2I1- 

The series (JTJ for I ^ then formally determine the lattice sums D^. Indeed, the free 
space Green's function is known to possess a partial wave expansion, 

G >,r,r') = -zACj^Jii^H+iar^yLir^yi^), (9) 

L 

where r> (r<) is the larger (smaller) of the |r| and |r'|. The numerical constant C [see 
Eq. ffTUj) or Eq. ()149|) below] is basically the prefactor in Eq. (J5J. Note that 

"-^■"-{jlE. <10) 

When the partial wave expansion is substituted for Gq in the series (jHJ) for G a, then, 
according to Eqs. (JJJ), (JHJ), one has 

D L (a, k) = -iCA l / 2 5 w - iAC ^ 'W/>r n )K(r n )e ik - r », (11) 

r„GA 
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where a prime on summation sign will here and hereafter indicate that the term r n = is 
omitted from the sum. Note that the series is, up to a constant term and a proportionality 
factor, the special case of series (JIJ for I ^ and R = 0. 

It is worthwhile to point out that one is often rather interested in the lattice sums Dl 
than at G a itself. For instance, since D^s do not depend on R, G a can be evaluated 
using Eqs. (J7J) and (JBJ) at any observation point with the same set of the lattice sums. 
This can in turn significantly speed-up numerical solutions of various boundary-value 
problems. Knowledge of the lattice sums Dl is a key to efficient numerical analysis 
of various ab-initio first-principle multiple-scattering problems, such as band structure 
calculation within Korringa-Kohn-Rostoker (KKR) theories 

1551 IdTH IHIl E21 EB1 EU ESj, and diffraction problems by periodic structures, or gratings, 
using the so-called layer KKR (LKKR) theories or equivalents thereof, either for quantum 
(electron) waves within low-energy electron diffraction (LEED) theory jlj |20J 122 EHU EI] 
or, for various classical (acoustic, elastic, electromagnetic, water) waves [5J EH1 EEJ 
E3 EH IZH H21 EH Ell- Lattice sums also arise in quantizing classically ergodic systems 
such as Sinai's billiard 24 and its various electromagnetic analogs [2U [75] . Moreover, 
the spectrum within the one-particle theory of periodic point (zero-range) interactions is 
determined as the set of those z which, for a given kii and a, satisfy an implicit equation 

nam 

Doo(iV=i,k||) = a. (12) 

The spectral parameter a here is a boundary-condition parameter which determines the 
asymptotic of eigenfunctions at the lattice points and is the same for all the eigenfunctions 

maim urn. 

3 Dual representations of quasi-periodic Goa 

Let A* be a corresponding dual (momentum) lattice, i.e., for any r n G A and k s G A* one 
has r n ■ k s = 27rA^, where is an integer. It is well known that Goa has an alternative 
representation as a sum over the dual lattice A*. For example, in the bulk case, i.e., when 
g?a = d, the dual sum representation of Goa is 

r i ^ nl 1 v- e" k+k - ) - R v- <fi(k + K, rfflk + k„ rQ 

0A< ' ' H) ^ St. °* ~ * + W = Si. °> ~ (" + W ' < 3> 




(14) 



normalized to unity in the fundamental (Wigner-Seitz) domain, 



I ip(k+ k^rWk + k^r) dr = 6 SS >, 
ws 



(15) 



Vq being the volume of a unit-cell of A. A dual representation is sometime called an 
eigenfunction expansion of Goa- Often the respective representations Q and (fT3j) of G a 
are also called the spatial- domain and spectral- domain forms, respectively [2E1 1213 EH1 E3 



In the following, the dual sum representation of G a in the quasi-periodic case of ID 
periodicity in 3D will be derived and its convergence properties will be discussed. However, 
before proceeding any further, it turns out expedient to provide some supplementary 
geometrical definitions which are up to minor variations adopted, for instance, within 
LEED and LKKR theories @1 12D1 1221 12EI IZ3 EH1 EZl EH1 EHl IZ01 CZ21 E31 EU - 

3.1 Supplementary geometrical definitions 

In the quasi-periodic case we define the respective parallel and perpendicular components 
rn and r_|_ of a given vector r = rn +r_|_ with respect to the <iA-dimensional plane containing 
the Bravais lattice A (i.e., line for d\ = 1 and surface for d^ = 2) and its normal n , 
respectively. Then for any r„ G A, 



The respective projections ky and kj_ of wave vector k = ky + kj_ are then defined 
in like manner with respect to A*. Obviously, in a quasi-periodic case the quantities 
G a(o~, k, R), D^(a, k, R) and D L (a, k) entering Eqs. (jSJ), (J7J), (jHI), (fTTj) are only functions 
of ky. Therefore, in the quasi-periodic case, i.e., if the underlying lattice A is of lower 
dimensionality than the embedding space (d\ < d), it is more appropriate to call merely 
the projection ky as the Bloch momentum. 

A plane wave e* k r incident on a scattering plane of identical scatterers arranged reg- 
ularly on A would (see Fig. P), in general, be diffracted (transmitted) to a wave with a 



35]). 



• r 



= 0. 



(16) 
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wave vector K n (K+), where = (ky + k n ,Kf n ), 

K± - ±K _ / ±^ " l k ll + knlT /2 , ^ 2 > Ik,| + K\ 2 (17) 

JX ±n ~ n=xv -Ln — 1 _l_„TlL-„ i 1, 12 ^211/2 _2 ^ |u„ , i, 12 I 1 '/ 1 



v ±z[|k||+k n | 2 -a 2 ] 1/2 , ° 2 < m+Ki ■ 
k n e A*, and |k| = |K„| = o. Here Kj_ n is indicated as a scalar, which is definitely 
true for d — d\ — 1. In the case of a ID periodic chain in 3D (d — d\ — 2) Kj_ n will be 
taken as the projection of Kj_ n on the plane spanned by the wave vectors of incident and 
diffracted beams. In the above definition, the projection Ky n = ky + k n is real but the 
normal projection K± n can be either real or imaginary. In the case of real K± n we speak 
of a propagating wave, and in the case of imaginary K± n of an evanescent wave. 

In the present case the scatterers are absent. Nevertheless, it turns out expedient to 
define wave vectors and the respective projections K||„ and K± n even in the free-space 
case. 

3.2 Resulting series 

In order to establish a dual representation of a quasi-periodic G a for d — d\ = 2, i.e. 
a ID periodicity along the x-axis in 3D, one first substitutes the integral representation 
(0) of Gq into defining equation J3J) of a free-space quasi-periodic Green's function G \. 
Then the Poisson formula 

_ e^~^ = ^ £ <f(qn - k,| - k n ) (18) 

is applied resulting in 

(n \d A -d r iq ± -Rx+i(k||+k n )-R|| 

GoA(ff , k|| , R) _ _J_ ^ j ________ dq± , (19) 

Now the 2D plane-wave expansion ()146|) is applied to e* q± Rx . Using the orthonormality 
of cylindrical harmonics YJ = e l1 ^/ \^2n one finds [Eq. ([145)1 ] 

Y l {<f>)d<t> = V2^6 l0 . (20) 



o 



Therefore, integration in the integral representation (fTHj) of Gqa over results in 

G A(a,k||,R) - 2_ e " / „2_„2 _| t — r 12 T ?V ^ 



V ° k^* ^ 



t~ _ e^ k li +k ")- R ii__- (^n|it! ± |). (21) 



4vo k„ eA * 



Here in going from the first to second equality the integral identity ()150|) for 2D case has 
been applied. Thereby the sum over r n G A has been transformed into a sum over k n G A* 
resulting in the so-called dual representation (spectral domain form) of G a- 



3.2.1 Complementary cases 

For completeness, in the case of codimension one (d — d\= 1), one first performs a partial 
integral over R^ in the integral representation (|SJ) of free-space Green's function. This 
amounts to picking up a residue of a contour integral in the complex plane according to 
Cauchy theorem resulting in 



m - r *qii-R|i+y^-q|||R±l 



Substituting the integral representation of Gq(ct, R) into defining equation (jSJ) of a free- 
space quasi-periodic Green's function G a then results in 



GV(cr,k||,R) 



iq r nu+iJa 2 -cifAR ± \ 

m / e " " v ii 



— T 



o"^ — q 

i(k,|+k„)-R||+ii<r_ Ln |i?_ L | 



k ||)- r ||n 

,r„£A 



dqil 



2t» ^ K 



J] j^^h^iK^R^). (23) 



2w 

where K± n is given by Eq. (j!7|) . Here in going from the first to second equality the Poisson 
formula (fTHj) has been applied. Note in passing that the respective dual representations for 
a ID periodicity in 2D and 2D periodicity in 3D are formally identical, the only difference 
being the dimensionality of A* in (|23J). 

Since exponentials in the second equality in (}2*3*j) can be rewritten as a product of two 
"eigenfunctions" , the dual representation can be recast in the form of an eigenfunction 
expansion of G a [cf. Eq. (jl3*))]. On the other hand, since H (K± n \R±\) in Eq. (j2~Tj) cannot 
be factored out as a product of two "eigenfuctions" , an eigenfunction interpretation of a 
dual representation is obscured for the case of a ID periodicity in 3D. 
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3.3 Convergence and limiting cases 

Note that, according to Eq. (|17jh K± n is purely imaginary with positive imaginary part 
for |k|| +k n | > a. For purely imaginary argument z = ix with x > the Hankel functions 
Hq (ix) are related to modified Bessel functions K (x) (Eq. (9.6.4) of Ref. [HI|), whereas 
(ix) in the series for d — = 2 are related to modified Bessel functions of third kind 
\A"/ (2x)Ki/ 2 (x) (Eq. (10.2.15) of Ref. [Ml). This results in rapidly decaying terms 
and exponential convergence (see Eqs. (9.7.2) and (10.2.17) of Ref. [HI]). Exponential 
convergence can also be directly inferred from the explicit expression for 

4 1} W = f- (24) 

One then easily finds that the respective terms in series (|23|) become exponentially de- 
creasing with increasing |kn + k n | for R± ^ 0. Since Hq(z) ~ \Jlj (ttz) e*^ 71 "/ 4 ) for 
\z\ — > oo and — n < argz < 2n (see Eq. (9.2.3) of Ref. [TH]; see also Eq. (J2J) above), 
similar applies to the series (|21|) . Therefore, although the convergence of a dual repre- 
sentation of Green's function is initially (for |k|| + k n | < a) slow as the series consists of 
mere oscillating terms, afterward (for |kii +k n | > a) convergence becomes exponential for 
R± 7^ (assuming as usual K± n ^ 0). 

Note in passing that all the dual representations of the reduced sums are nonanalytic 
in R± as they are functions of \R±\. In the case of series one has 

\R±\h { o\K ±n \R ± \) -> -i/K ±n as \R ± \ -> 0. (25) 

Therefore, absolute convergence of the resulting series in the limiting case cannot be 
established even for a ID periodicity in 2D (jH], PP- 51-52). Even worse, in the case of a 
ID periodicity in 3D individual terms of the series (J21|) possess a logarithmic singularity 
(see Eqs. (9.1.3), (9.1.13) of Ref. HEJ). 

To this end, it has been demonstrated that absolute convergence of dual representa- 
tions can only be established under the assumption of R± ^ 0. Additionally, obtaining 
the quasi-periodic Green's function of the Laplace equation from that of the Helmholtz 
equation by taking the limit o — » in the resulting expressions (see Sec. El below) is 
problematic when starting from a dual representation. Since 

H+{K ±n \R ± \) ^H+{*\k\\+ k n\\R±\) asa-^0, (26) 

11 



absolute convergence can again be established only for R± ^ 0. In order to resolve the 
above problems it turns out expedient to invoke representations which converge uniformly 
and absolutely with respect to R. Such representations are known as the Ewald represen- 
tations [T%l ll9|. Additional bonus of the Ewald representations is that they enable one to 
investigate analytic properties of quasi-periodic Green's functions in the complex variable 
z = a 2 . 

As a final remark of this section note that dual representations in the quasi-periodic 
case can also be established by applying an Ewald integral representation of Green's 
function and the generalized Jacobi identity. This path, which has been originally followed 
by McRae j3H] for a 2D periodicity in 3D, is outlined in Appendix lEl 

4 Ewald representations of quasi-periodic Goa 

In going from the spatial domain form [Eq. ©] to the respective spectral domain forms 
of Goa [Eqs. (}2~T]) . (J23j) ]. the summation over A has been fully replaced by a summation 
over A*. In this section, starting from the spectral domain forms of Goa a half-step 
backward will be performed resulting in a hybrid Ewald representation of Goa- The 
Ewald representation of Goa involves sums over both A and A* and, in contrast to a dual, 
spectral domain, form of Goa, is valid for all R± and uniformly convergent with respect 
to bounded sets of R, provided that R ^ A. 

In deriving the Ewald representation of Goa, we first recall formulae (10.1.1) and (9.1.6) 
of Ref. [12] and recast Hq in the series (|2^|) as 

4 1 'w = (^) 1/2 <w = -i© 1/2 «%w- < 27 > 

Therefore, a dual representation of G a in any quasi-periodic case [Eqs. (|2*Tj). (j2Hj)] can 
be rewritten as a sum of cylindrical Hankel functions of an appropriate order. 
Now for d — d\ = 2 we shall introduce 

Gj,a(c, k, R) =~4- (^y /2 H%(KUR ± \)e^^. (28) 

In the following, G^a will be called an analytic form of Goa- Obviously Goa = Gva\ u= q- 
Next it turns out expedient to employ the following integral representation [see Eq. ()138j) 
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of Appendix |Xj , 



oo exp iq 



H%(K ±n \R ± \) = ± " C /^ e H K ^-\n,\VO dC (29) 

K±n J ~ 7™ Jo+ 



The lower limit + indicates that the contour integral starts from in the direction of the 
positive real axis. Here we have used the convention (see Ref. jlH], p. 589) that 

I D i \ *V2 ( 
it i \ \ v 



K^J 6XP 2 ,ln 



R\ 



iargK ±n } , (30) 



where the argument of K± n takes the values or n/2, and cj) n is given by 

(j) n = ix - 2aigK ±n . (31) 

From Eq. (J2*9j) it follows that Hq (K± n \Rx\), and hence also Goa ; can be analytically 
continued in the complex parameter v. This explains the reason why G u \ defined by Eq. 
PHj) has been called an analytic form of Gq\. Indeed, as soon as Re v > 0, 

H*il /2 (z) ~ -i——T{v/2){z/2Y u l 2 as z -> (32) 

(Eqs. (9.1.6) and (9.1.9) of Ref. [16J). Therefore, all the terms in the series ()28|) for 
Re v > are singularity free as — > 0, and the original logarithmic singularity of 
Hq (K± n \R±\) in the limit \R±\ in the dual representation (|2*T|) of Goa is thereby 
avoided. Additionally, the series in (|28J) can be easily seen as an analytic function of v 
for all values of R, provided that Re v > 2d\. In the latter case, upon using asymptotic 
(|52~|l . the elementary products \R±\ v ^ 2 H_ v , 2 {K±_ n \R±\) in the series ()28|1 can be uniformly 
bounded for all n by a finite number as — > 0. On the other hand, the asymptotics of 
H { }'l /2 (z) = e™/ 2 H^ 2 (z) (Eq. (9.1.6) of Ref. as K ±n — > oo is determined according; 
to Eq. (0) with d = 2. Consequently, the series ()28|) can be uniformly bounded by the 



series 

K ±n 2 for a11 R ' N ° W tlle 

series K ± ^ 2 is absolutely convergent for Re v > 2d 
(Ref. [IB], pp. 51-52). Therefore, since the sum in Eq. (|2*B*|) is absolutely and uniformly 
convergent for all R, it defines an analytic function of u for all values of R± and Rii if Re 
v > 2d A (Ref. HE], Sec. 5.32). 

Now, the task is to find an analytic continuation G vA from a domain Re v > 2d\ to a 
domain containing v — 0. Here it has been implicitly assumed (and will be proved later 



A 



13 



on) that the analytical of G v \ defined for Re v > 2d\ by the series ()28|) will yield our Goa 
defined by Eq. Q. The necessity of an analytic continuation in the quasi-periodic case 
makes derivation of the Ewald representation of G*oa fundamentally different from that in 
the bulk case. This analytical continuation procedure is analogous to finding the values 
of the Riemann ^-function outside the domain of absolute convergence, Re v < 1, of its 
defining series 

oo 

C{v) = Y. n ~ r ( 33 ) 

n=l 

In order to find out an analytic continuation G v \ in a domain containing v = 0, one 
substitutes back into (J28|) which results in 



1 r—oo 

G„ A (<7,k|.,R) = — y e l ( k ii +k "> R ii / C /2~i e h(Kl n C-\R^/O dC 

4 ™°kfil* V 

— - V e i(k ii +k " ) - R ii / C w/2-1 e5 (^InC-lflxlVO d£ 
4 ™° kfe. ^ 



(34) 



The first sum has only a limited number of terms so that the order of integration and 
summation can be inverted. By a straightforward generalization of Riemann's method 
(see Ref. [IB], p. 273), the inversion of integration and summation also holds for the 
second sum in (|34|). provided that Re v > 2d^ and if always Re ( > on the contour of 
integration [2T] . 

The remaining two steps in the derivation of the Ewald representation are essentially 
those used by Epstein in an analytic continuation of his zeta functions fUl : 

• the resulting contour integral is split in two parts by taking a point rj somewhere in 
the domain Re r] > 0, \r)\ < oo. 



the generalized Jacobi identity ([163)1 . which is valid for Re £ > 0, is applied for the 
part of the integral from to rj with ( = l/(2£ 2 ), and r s = — r s , yielding 

E -(k n +k||) 2 C/2+i(k||+k„)-R|| _ ^0 „-(R||+r„) 2 /(2C)-ik r r n /oe\ 

kn£A* v SJ r„eA 
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Following the first step, G u \ is expressed as the sum of two terms, 

G vA (a, k,|, R) = G 1 (a, k„, R) + G 2 (a, k,,, R), (36) 

where the respective G\ and G2 contributions result from the respective contour integrals 
over (0,rj) and (77,00). Obviously, although each of the partial integrals depends on 77, 
called the Ewald parameter, their sum does not. 

Following the second step, Eq. (J31j) is transformed into 



-1 p— OO J-n' 

GuA CR) = L_ / y ei (k ||+ k„).R |ir /2- le i(Ki nC -|iixl 2 /C) dC 

/ e i(k||+k n ).R|| (:I ,/2-l e I(iCi n C-|«x| 2 /C)^ 

4tTW Jr, ^ 



(2^3/^2 j y 2-i c c ' 



(37) 

r„eA 

The latter expression is an analytic function of v for all values of v if |R| 7^ 0, or more 
generally if R ^ A, and it represents the sought analytic continuation of (|28|) for Re 
v < 2d\. Note in passing that for Re A analyticity can only be established if Re v > 1. 
Otherwise, the last integral diverges for r n = — R. 

On putting v = 0, and hence assuming R ^ A, inverting again the order of summation 
and integration (since it can be allowed), and, substituting ( — > l/( in the last integral, 

"OO exp «</ 



1 /-ooexpztp,! 

GW^WR) = — Y e^ii+^ii / ' C -i e K*inC-l^l70 dC 

- ^ f-) V2 F ^ ik|| - rs /°° C - 1/2 eM- 2 /C-(K + r s ) 2 c] dC (38) 

The restriction Re ( > going back to the integral representation ()29jl can be now 
removed. 

4.1 Complementary cases 

For completeness, for d — d\ = 1 one would begin with the analytic form 

GM,R) = (f)" 2 E, (f^ ^» /2 (A-,„|fixl)e*.*>-, (39) 
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Similarly as in the preceding case, G v ^ defines for Re v > 2d\ an analytic function of v 
for all values of R± and Rh, provided that it remains R ^ A. The task is now to find an 
analytic continuation G v \ from a domain Re v > 2d\ to a domain containing v — 1. After 
repeating the steps which led from Eq. (J25|) to Eq. (JHHj) . on putting v = 1 and inverting 
again the order of summation and integration, one finds for d\ = 2, d = 3 jlj |2*U] 

1 1 /2 Z* 00 CX P ^ n 

CoA(a,k||,R) = —(-) T e J ( k ii +k "> R " / r i/ 2e |(xi„ C -|«xlVC) rfC 

"A (-) 1/2 E e " k "' rS r ^ 1/2e " [<7VC - (R+rs)2c] dC- (40) 
47r V2/ _/]_/„ 



r s £A 



Similarly, for d\ — 1, d — 2 one arrives (see, e.g., Ref. (2H|) at 

1 , _ x 1 /O Z" 00 ex P *</»n 

G 0A (a, k|| ,R) = -^-(|) ^e^ +k ")- R ii/ C^^LC-I^VC) dC 

1 /*oo 

__L Vy*""'' / r'eH^/C-CR+r^c] dC (41) 

It can be proved directly (see, e.g., Appendix 3 of Ref. [21]) that Ewald representa- 
tions (jHHI); ()40|) . (|4*T|) satisfy Eq. (£Q) and the boundary conditions. Therefore, they are 
required Green's function, which provides a posteriori justification of the outlined analytic 
continuation procedure. 

Since the respective dual representations for a ID periodicity in 2D and 2D periodicity 
in 3D are formally identical, the terms involving a sum over reciprocal lattice in Eqs. ([4*U|) 
and ()41|) are identical. Surprisingly enough, the terms involving a sum over direct lattice 
in Eqs. ([3~%j) and (j4*nj) are identical too. 

4.2 Ewald vs dual representations 

One has (see Appendix 1 of Ref. [2*T] ) 

f»oo exp icf> n 

£-i/2 e i(^i n C-|^l 2 /C)^ = 

, e iK ±n \R x \ roc 



and 



tK ±n \H±\ r<x> 

-V^F 6 . - [ r 3/ M (^AH^I 2 d( 
iK±n Ji/n 

poo exp i(p n 
Jr) 
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poo 

7riH^(K ±n \R ± \)- / c^^'i^-^Uc. 

Jl/n 



'1/7, 

For K\ n > this can be shown by deforming integration contour in Eqs. ()38)1 . (}4T)|) . fiT| 
to that shown in Fig. |2k and upon invoking Jordan's lemma for the integration along 
quarter circles (|48j. p. 115). For Kj_ n < one then takes the contour as shown in Fig. 
12)3. Therefore, a comparison of dual representations (j2Tj) . (j2*Hj) of quasi-periodic free-space 
Green's function with their respective Ewald representations (j2EJ), (j4Ti|) . shows that 
to each term of a dual representation the second term (integral above) is added to make 
the series convergent uniformly with respect to R± (provided that R ^ A). These terms 
are then compensated by the series over A. 

For a sufficiently large rj, G a can often be well approximated by the series over A. 
This approximation to Goa is called the incomplete Ewald summation |54j . 

5 Calculation of the lattice sums Di 

The lattice sums D L have been defined by Eq. (jSJ) as the expansion coefficients of G a 
in terms of the regular (cylindrical in 2D, spherical in 3D) waves, or, alternatively, as 
the Schlomilch series (|llj). Analytic closed expression of the lattice sums can only be 
obtained in the particular case of / = and a ID lattice A with a period a in 3D |T4~1 115j. 
Indeed, assuming the elementary identity 

ln(l -*)=-£-, (42) 



one obtains 



, n 

n>0 



k r n a \n\ 

r„6A n^O 1 1 

= _I {In [1 - e ia{ ^0] + In [1 - e -(-+ k ll)] } 

= -- In \e 2iaa - 2 cos(k|,a)e ia<T + ll . (43) 
a 

Hence, upon using that — iACY 00 = —ia/^/An in 3D, 

Ax>(<7,k[i) = -i-^ + ^^ln[e 2ia<T -2cos(k||a)e w,T + ll 
V47r v47ra 



! In {2 [cos(cra) - cos(k||a)] } , (44) 



Ana 
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which is up to the prefactor of — l/(\/47ra) the 7-function of Karpeshina (see Eq. (29) 
of Ref. [T3| for a = is). (When solving Eq. (|12|). the principal branch of logarithm 
is assumed in the above expression for Dqq and our spectral parameter a has also been 
rescaled compared to that of Karpeshina [T3"l l4T] with the above prefactor.) In accordance 
with our notation, boldface r n and ky are numbers which can be either positive or negative, 
whereas r n > stands for absolute value. (It is reminded here that the energy operator in 
the one-particle theory of periodic point interactions is constructed in terms of the operator 
of multiplication by -Doo (7 function of Karpeshina JH1E1) an d that Dqo determines the 
spectrum according to Eq. (|12p.) 

Invoking that Dl's are independent of R, the lattice sums in all remaining cases are 
calculated as [H 1201 

D L (a,k l{ ) = lim v— \— (fyi(k)D A (a,k h R)d^ (45) 

where § dVL-& denotes the angular integration over all directions of R. In calculating 
Dl, the respective Ewald representations (jHEJ), (f4*Tf of Goa are substituted in the 

defining Eq. (JZj) for D^(a, ky, R). The two series in the respective Ewald representations 
(138)) . ()4(J|) . (|41|) are uniformly convergent with respect to R so that the series can be term- 
wise integrated when is calculated according to Eq. ()45|) . Following a hybrid character 
of the Ewald representations (jSHJ), PHJ) . (jUJ), the respective Dl are conventionally written 
as a sum H EI3 EEl EH 

D L (a,k {l ) = D^ia,^) + D^(a,k\\) + D L 3 \a), (46) 

where (D^) involves a sum over reciprocal lattice (all r n 7^ terms of the direct 

lattice). D L 3 ' is the term which combines Gq(R) and the r n = contribution of the direct 

(3) 

lattice sum G2. D L is only nonzero for 1 = 0, 

Df = D®6 L0 . (47) 

In the following, the respective contributions D L l \ D L 2 \ and D L 3 ^ will be calculated. 
For reader not interested in an explicit derivation of results, the resulting expressions are 
given by Eqs. flTU^J) . and IfTTBjl [see Eqs. JTU1, dSJ for d A = 2, d = 3 and 

Eqs. dEHD , (UnnD, (gH for d A = l, d = 2]. 



18 



5.1 Consequences of the reflection symmetry for the lattice sums 
in the quasi-periodic case 

Assuming standard spherical coordinates, one has 

F L (R|I - Rl) = (-l)' +m F L (R|| + Rj_). (48) 

Therefore, for the lattice plane perpendicular to the £-axis 

D L = 0, l+m odd (49) 

for both c?a = 1) d — 3 and = 2, d = 3 cases. This identity follows upon combining the 
property (|184j) with the expansion (jSJ). In fact [see Eq. (JTTj) and Eqs. ()67|h (jl(JU|) below], 
it will be shown that (for the lattice plane perpendicular to the z-axis) the property (J49j) 
holds for each of the contributions j = 1, 2, 3, separately. 

5.2 Calculation of 
5.2.1 General part 

As it has been alluded to above, the contribution Dj*/ derives from the sum over reciprocal 
lattice A in the corresponding Ewald representation of Gq\. According to the Ewald 
representations (|38|). (|4U|h (jTTj) of the quasi-periodic Green's functions, 

DV = lr^%\ (50) 

L 2w (2tt) c l ' v ; 

where 

j(D = l im _J_ I y L *(R)e J ( k ll +k ^- R H / r c e I(^ sC -|fixl 2 /C) dCdQ ^ (51) 

R -*° j\i\{^R) J J v 

and l/2<c=(d — d^)/2 < 1. The exponential decrease of the integrand with increasing 
C for the integration over ( (assuming as usual Kj_ s ^ 0) guarantees that the order of 
integration can be inverted. In order to perform the latter integral, the exponential is 
expanded into a power series of | R± | 2 resulting in 



lim — 



x ( lim 1 - & F^(R)e^ k n +ks > R n|ii: ± | 2n dn R 
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n=0 

2l ^ (1 " cW E||r(i- c - n , e -^^M\ K ^-i\ (52) 

n=0 ' ^ ' 



where 

Jj(n) = lira l — ^— /F;(R)e^i +k ^' R ii|i?J 2 "^ R , (53) 

and r is the incomplete gamma function (see Eq. (6.5.3) of Ref. |16|). In the second 
equality in (J52|) we have used in the integral over £ the substitution 



2p' 

C = ^5-*» (54) 



s 



which leads to 



X i'Xp ;r,)„ / 2 ^ 1~ C~ ™ 



sJ J e-™Kl s r)/2 

1—c—n / T r2 



In the final equality in (J52*|l we have substituted 



(55) 



K 2 \ n+c-1 f_-\\n 

e -*i±Ls) = gd-c^i-cLiL^ 1 ). (56) 

2 / 2 n ±s v ; 



Now upon combining Eqs. (JH0|) and ()52jl 



U k s 6A* n=0 



x ( lim } </y L *(R)e i(k ii +kB) - R ii| J Rj.| 2n df2 R ) . (57) 

It order to finish the calculation of it remains to perform the angular integration 
and the limit R — > 0. In the following this limit will be provided for various particular 
cases. 

5.2.2 The case of a ID periodicity in 3D 

If the ID lattice is oriented along the x-axis, Rj_ = Rcos6*, Ry = Rsin#cos0, and hence 
| HL_i_ [ = R\ cos0|, (k|| + k s ) • R|| = |k|| + ksl^Rsin^ cos ^0k||+k s — 0^ , (58) 
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where <p u is the polar angle of the vector u in the plane containing A* (A). The only 
difference with respect to the case of a 2D periodicity in 3D, which has been treated by 
Kambe [3J ETJ E] is merely in that the values of 0k||+k 3 are no longer from the interval 
[0, 27r) but are restricted to either or n. According to Eq. (J32J), one has 



(59) 



r(l) _ ,-m-M A7 ,f 1 F ( -iri K lsV \ K 2n ( y 

h ~ l ^^WH ' 2 ) K±s {R^oj l (aR)J 
where Iq involves the following angular integral (dfl = sin OdOdcf)) , 

ig = f\in6d6P\ m \ C os6)( C os6) 2n f 2 \-^\^s\Rs i necos(^ +ks - (t> ) # (6Q) 
Jo Jo 

Integrating first by <p one finds, 

/*7T 

= 27ri H e"^ k n +ks / P, H (cos#)(cos#) 2n J H (|k|| +k 3 \RsmO) sin0<J0, (61) 
where we have used that 

e -*"*W*«*(*>-0 c/0 = 2m |m| e- m * Jh(^). (62) 
The latter identity can be derived from (cf. Eqs. (9.1.44-45) of Ref. [TT)] ) 

oo 

e i2COS *= ^ zl'lj| Z |(z)e^. (63) 

Since J„(z) is an entire function of z, the Bessel function in Eq. f|61j) can be expanded 
into power series of its argument (see Eq. (9.1.10) of Ref. [TTj]). 

Afterward Eq. (joTj) becomes one obtains 

jn _ o„>| -«^ |+ k s yV , y 7 - [l k [| +k,li?]H+2j 



2M+2jj!(| m | + j)l 



(65) 



where 



4 = J (cosfl) 2n (sinfl)H+2i+i p\ m \cos9)d6 = j x 2n {l - x 2 ) J ' +H/2 P, H (x)cfe. (66) 
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It will turn out that, in addition to the general property flOj) . each term Z?j*2 = 0, 
j = 1, 2, 3, for I — \m\ odd. In particular, one can show that the integral ()66j) vanishes 
unless I — \m\ is even, i.e., 

D im =0, I- \m\ odd. (67) 

Indeed, 

pH (cos = (-lj'-HpHjcosf, - 9)]. (68) 
Now, if F is a real function such that F(8) = F(tt — 6), then for / — |m| even one finds 

F(6)Pi ml (cos 0) d# = 2 /" /2 F(0)P' m| (cos 0) dO. (69) 

On the other hand for / — \m\ odd the integral vanishes. Since in our case 

F(9) = (cos#) 2n (sin#) H+2j+1 = F(tt - 6), 

it follows that F e = for i — \m\ odd. 
For / — \m\ even 

l im J^£L ^o, j+n> . (ro) 

Therefore, one only needs to investigate the case 

I — Iml , , 

3+n<^-^- (71) 

Upon expanding (1 — x 2 ) J in (jo^j) according to binomial theorem, P e is rewritten as 

5 



lJ o = ( 1) f - x 2 ) ]ml/2 Pi ml (x)dx. (72) 

According to Eq. (2.17.2.7) of Ref. [77j, the integral 

J 1 x t {l-x 2 )\ m \l 2 P\ m \x)dx (73) 

vanishes unless 

\m\ < I < t + \m\, (74) 
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or, in our case, unless 

I — \m\ 

\m\ < I < 2n + 2s + \m\ n + s > — - — • (75) 

Upon combining the conditions (j7ip and ()75|) . it follows that the only nonzero contribution 
to I 3 g in the R — > limit arises when simultaneously 

I — \m\ 

s = j and n + j = — - — ■ (76) 
In the latter case, Eq. (2.17.2.6) of Ref. [77j implies 



2 m (2n + 2j)!(/ + 


m\)\[n + j + (I + \ 


m\)/2\\ 


[n + j - (/- 


m|)/2]!(2n + 2j + 1 + 


m\ 


+m-\ 


m\)\ 



= (-iy^g±Ml. (77) 

1 ; (2/+1)! 1 ; 

Here in the last equation we have substituted for n + j according to Eq. (JTHJ). (Note 
in passing that Eq. (2.17.2.6) of Ref. [77] differs by a factor (— l) m /2 compared to Eq. 
(7.132.5) of Ref. [ZHj due to a slightly different definition of associated Legendre functions.) 
The constraints ()76|) imply 

\m\ + 2] = 1- 2n, (78) 

and 

n < L_L ■ (79) 



Therefore, the sum over n in I L becomes a finite sum. 



Consequently, as R — > 0, 



H l n ~ e ii (-1J J 



2i-an[(/_ | m |)/2-n]![(Z + |m|)/2-n]! 



, . 2^/!(/ + |m|)! . 
^ J (2/ + 1)! 



2'" 2 "[(/- |m|)/2-n]![(Z+ |m|)/2-n]! (2/+1)! ' l ; 
Now (Eq. (10.1.2) of Ref. PU) 

J,MH (2iTi^ M) (jR ^ 0) - (81) 
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Therefore, in the limit R — > 0: 



r(l) 



V ■ J, 

2v^F— [(2/ + \m\)\(l - |m|)!] 1/2 e~ mu ^^ 



(/-|m|)/2 

E r 



n=0 



-n, e 



\k l{ +K\/cT] l - 2n [K ±s /a 



12n 



n![(Z - |m|)/2 - n]![(Z + |m|)/2 - n}\ 



12) 



Consequently, for I — |m| even, 



1 Z m 



— =-^[(21 + \m\)\(l - \m\)\f 2 e 

2 V^o 2' ^ 



n+k s 



(J-M)/2 

E r 



n=0 



-n, e 



|k||+k s |/a]'- 2 "[K ±s /a] 



I2n 



ra![(Z - |m|)/2 - n]![(Z + |m|)/2 - n] 



whereas [Eq. (JB7J) ] 



= 



/ — I m I odd. 



(83) 
(84) 



5.2.3 Complementary cases 



For d\ = 2 and d = 3, provided that lattice plane is perpendicular to the z-axis, one can 
repeat most of the steps presented here, the only change being c = 1/2 instead of c = 1, 
and thereby arriving at the Kambe's expression [201 1221 EZ] 



£>LV,k|| 



rW 



2w V27r 

1 r m+1 



ctoo 2' 
(i-IH)/a 

E r (V2 



[(2i + l)(( + |m|)!((- |m|)!] 1/a yj e"" 1 1 



n=0 



n, e 



1/2 „" im ^ 
k s 6A* 

[Ikii+k^l/a]'- 2 "^/^ 2 "- 1 
n![(Z - |m|)/2 - n]![(Z + \m\)/2 - n}\ 



15) 



For d\ = 1 and d = 2 one then finds j2H] 

[|»l/2] 



A (1) (a,k|| 



y/2(JV 



W v 'f j_ r ( 1/2 _ n e -*?M) 



k s eA* n=0 

|k|| + k s \/a]^~ 2n [K± s /a]^^ I exp 

exp 



(|Z|-2n)! 



2 / 

-*(|/|-2n)0 k|| + k . . />() 
- 2n)0 k||+ks . / < 0. 



(86) 
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where vq is now the length of the primitive cell of A and [\l\/2] stands for the integral part 
of \l\/2. 

5.2.4 Convergence and recurrence relations 

In virtue of the asymptotic behavior 

T(a, z) ~ z a ~ x e~ z as z -> oo, |arg z\ < 3vr/2 (87) 

(Eq. (6.5.32) of Ref. ^H]) it is straightforward to verify that convergence of the series 
on the r.h.s. of expressions (JBHJ), ([55]). (jHEJ) is exponential for sufficiently large K± s . One 
can verify that the lattice sums are dimensionless for d = 2 whereas for d = 3 the 
lattice sums have dimension [1/length] [see discussion below Eq. (|191|) ]. From the 
computational point of view, the incomplete gamma function in final expressions 
(|85j) , dHBj) can be derived successively by the recurrence formula (201 1 

bT(b,x) = T(b + l,x) -x b e~ x 

from the value for n = 0: 

r(l/2,x) = <{ v / ^-2/ ^e-^t = v yerfc( v ^), argx = (gg) 

V / 7r±2i J Q e* dt, argx = =F7r. 



5.3.1 General part 



5.3 Calculation of D L 



As it has been alluded to above, the contribution D L derives from the sum over direct 
lattice A in the Ewald representation of Gqa, but with the term r n = excluded. According 
to Eq. (|4*5Jl taken in combination with the representations (f3~8]). (pfUj) . (}4*T]) of the quasi- 
periodic Green's functions, 

Df = - ( ^-lf, (90) 



where 



Here c = 1/2 for d = 3 and c = 1 for d = 2. Using the plane- wave expansion ([146)1 . which 
is also valid for complex arguments (see, e.g., Appendix 1 of Ref. jl] or Appendix A of 



25 



review by Tong [67)]). 



-(R+r s )^/2 = e -(R 2 +rlK/2 e i(i(r s y-R _ 

L 



Ae -(R^ )( /2 «r s i?) y* L (i s )y L (R). (92) 



The absolute value |/| here is only relevant for d = 2 and provided that the range of 
angular momenta is taken to be — oo < / < oo. 

Substituting (|H2^ back into equation for I L , taking the limit pHjl. and using that in 
any dimension 

Ji(az) a 1 

as z — > 0, (93) 



The substitution 



Ji/ V 



(94) 



(=—, r s 2 C = ^, dC = — a=^- (95) 

s 2u sS 2m ' s 2u 2 ' 2 K ' 

then yields 

jf = A(-lf2-^- l+c a 2 - 2c (ar s fy* L (r s ) [° u c ~^~ 2 e u ^ r '^ du. (96) 

To this end, the formula for Dj/ is valid in any quasi-periodic case. 
5.3.2 Quasi-periodic cases in 3D 

In 3D one has c = 1/2 for both ID and 2D periodicities. Eq. ()96|) then yields 

if = 47r(-iy2- / - 1 /V(ar s )'yi:(rg / V^^M/W du . (97) 

Jo 

Now, using representations (|38p. ()4L))1. 



r s eA 

where prime in ^ ' indicates as usual that the term with r s = is omitted. Since, for 
I + m odd, 

Y lm {ir/2, <f>) = 0, (99) 
Df = 0, l + m odd. (100) 
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(We recall here that the periodicity direction(s) has been assumed to be perpendicular to 
the z-axis.) This confirms [see Eq. (|67jl ] that, in addition to the general property (|49|) . 
each term Dj^ = 0, j = 1,2,3, for I — \m\ odd. In the remaining cases for I + m even 

M2 6) - ( l)WH)/» ijg [(2Z + l)(Z-|m|)!(f+H)!]^ , 

y^7r/2,0)-( i) [a-H)/2l!fa + M)/2l! e • (1U1) 



Therefore, 



imV ' 117 " 4tt 2 2 < 



[(2/ + l)(i- 


m|)!(/ + 


ml)!] 1 / 2 


[('- 


ro|)/2]![(Z + 


m\)/2\\ 



x V ' e - 4k i|- rs - im ^(ar s )' / w-'-^exp 



cr 2 r 2 

u - 



4m 



du, (102) 

/ . . Ml 

r s eA 

where a = a 2 rj/2 [see Eq. ()95|)]. Hence, the respective expressions for for g?a = 1, 
d = 3 and for d^ = 2, d = 3 1221 12Z| are formally identical (provided that periodicity 
direction(s) is (are) perpendicular to the z-axis), the only difference being in the lattice 
dimension. 

5.3.3 Complementary case of a ID periodicity in 2D 

In the remaining case for d\ — 1 and d = 2, in which case c = 1, one finds by a slight 
modification of the preceding derivation |2Sj 

£>{ a V,k|.) = - , ( r 1} !^ T'e- ik ^(ar s fe- il ^ 
1 V V 2l'l +1 v/2¥ 

v r s GA 



X 



^° M-I'l" 1 exp (u - ^pj du. (103) 



5.3.4 Convergence and recurrence relations 

Similarly as in the case of \ convergence of the series for on the r.h.s. of expressions 
fll()2j) . ()103j) is exponential for sufficiently large r s . Indeed, integrals 

U\i\ = [ u H'|-2 e «-* 2 r s 2 /(4u) du ( 1Q4 ) 
JO 

are finite integrals. Now, for \crr s \ > \l\ + 2 — c, the integrand is monotonically increasing 
from zero to ^a 2 r]/2) c ^^^ 2 e a2ri ^ 2 ~ r ^ ( - 2v ^ on the integration interval. Therefore 

|(<7r s )l%| < a 2 ^- 1 ^\r ] /2) c ^~ 1 r^e a2v/2 - r2j{2v) . (105) 
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From the computational point of view, the integral on the r.h.s. of the resulting 
expressions can easily be performed by a simple recurrence. Indeed, knowing the values 
of Uq and U\, the respective integrals U\i\ can be determined using recursion relation 

^Y) 2 U\l\ + l = + 1 - C)U W - U\ lhl + a -W-14* e «-« a rJ/(4a) > (1Q6) 

The recurrence here follows, as suggested by Kambe (see Appendix 2 of Ref. jl]), from a 
simple integration by parts. 

(2) 

As a consistency check, note that for d = 2 the lattice sums D L are dimensionless, 
whereas for d = 3 the lattice sums .D^ have dimension [1/length] [see discussion below 
Eq. (HH}]. 



5.4 Calculation of Df 
5.4.1 General part 



According to Eq. (|47jh the only nonzero term is -Dqq , or, for the sake of notation, D l ^ > . 
In 3D, for both ID and 2D periodicities, the term is calculated as the limit 



,(3) 



D ( q ] = Inn 



Similarly, in 2D case 



1 /7T\V2 



4tt 2 V2 



r i/2 e i(- 2 /C-R 2 C) rf( - Gl(a, R) 



(107) 



L>J 3) = lim 



4tt A/*, 



(108) 



It is reminded here that G^(a, R) in Eqs. (|1U7|) and ()1U8|) is the corresponding singular, 
or principal- value, part of Gq (a, R) [see Eq. © above]. 

In any dimension 3^o — Ar 1 ! 2 [Eq. ()144|) ]. where A is given by Eq. ()147j) . and (see 
Eqs. (9.1.12), (10.1.11) of Ref. |E|) 



J (<tR)->1 (R^O). 



(109) 



(3) 

Hence, calculation of Dq requires to perform integral 



roc 

j(3) = / C c e \ ( CT VC-R 2 C) d( 
Ji/n 



;no) 



28 



for either c = 1/2 (d — 3) or c = 1 (d — 2). Expanding e " 2//( - 2 ^- ) into power series and 
integrating term by term yields 

°° „2n poo 

lW = T,h C'-e-WdC (111) 



n=0 



Since all terms in the sum are positive, the exchange of integration and summation is 
justified whenever one shows that one of the sides exists (this will be shown later on). 
Upon substituting ( = 2t/R 2 , 



7 (3, = (flV2r .£l£*L / t -o-n e -, dt 



n=0 

where as usual T(a,x) is an incomplete gamma function (see Eq. (6.5.3) of Ref. [TB]). 
According to Eqs. (jl(J7|) and (|1(J8|) . we are only interested in the limit R — > 0. For 
n > 1 — c, which given c = 1/2, 1 translates into n > 1, one can integrate by parts, 
yielding 

/■oo i / p2 \ -c-n+1 

Av(2,) n + c-l\2 v J 1 J 1 ; 

Therefore, 

7 (3) = y — i — ^! + lim(jR 2 /2)c -i r r c e -t dt (113) 

A further treatment differs in different dimensions. 
5.4.2 Quasi-periodic cases in 3D 

In 3D one has c = 1/2 for both ID and 2D periodicities. Hence, 

oo / /-oo f R 2 /(2 V )\ 

R 2 /(2»?) JO J 

= T(l/2)-2-^ TJ - 2 e- R2 ^ + 0(R% (114) 
where T(l/2) = yfir. This asymptotic is consistent with that obtained by writing 

(115) 



T[l/2,R 2 /(2r])} = [ t- 1/2 e~ l dt = v^erfc (R/y/2rj 

Jn 2 /(2n) ^ 



'R 2 /(2 V ) 
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(see Eq. (6.5.17) of Ref. [TE]) and using that 



erfc z = 1 - -J= + 0{z 2 ) as z -> (116) 



/7T 

(see Eq. (7.2.4) of Ref. US]). Therefore, 

(j A ( ( T 2 )? /2) n - 1 / 2 



/(3) = _^ + " y IgJ^ _ +o(#) asi^O. (117) 

n!(n- 1/2) 



ra=0 



Collecting everything together back to Eq. (|l()7|l . the first term in J® cancels against 
Go(cr, R) leaving behind 

(U8) 

It is emphasized here that the respective contributions D lr ^ for a ID periodicity in 3D 
and a 2D periodicity in 3D (see Refs. [20112211211) are identical. 

5.4.3 Complementary case of ID periodicity in 2D 

For d\ = 1 and d — 2, in which case c = 1, one finds [2H] 



A (3) (a) = 1 



2x/2% 



7 + ln <.V2) + £ W 



n!n 

n=l 



<f ro = -4=Ei(<7V2)<fo, (H9) 

ZV Z7T 



where Ei is an exponential integral and 7 ~ 0.57721 56649 is the Euler constant (see Eqs. 
(5.1.10) and (6.1.3) of Ref. ^Ej, respectively). Note that, for A oriented along the x-axis, 
Di = Dzl, j = 1, 2, and hence, D\ = D_i, in accord with the fact that Goa only depends 
on y via \y\ [33"1 135]. 

5.4.4 Convergence 

Unlike preceding cases, convergence of series on the r.h.s. of Eqs. ()118|) and ()119|) is 
even faster than exponential. This can easily be verified by using Stirling's formula (Eq. 
(6.1.37) of Ref. [E]) 



(n+ 1)! ~ v / 2^n n e" n as n -> 00. (120) 

Again, as a consistency check, note that for d = 2 the lattice sums Z)^ are dimensionless, 
whereas for d = 3 the lattice sums D L have dimension [1/length] [see discussion below 
Eq. ()191j) ]. It is reminded here that rj in the above formulae has dimension of [length] 2 . 
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6 Laplace equation 



The quasi-periodic solutions of the Laplace equation are used to describe potential flows 
in fluid dynamics between parallel planes and in rectangular channels jHOl- Indeed, a 
Green's function representing a point source and satisfying the respective von Neumann 
and Dirichlet boundary conditions on a flow channel walls can be written as a sum and dif- 
ference of an appropriate Goa (corresponding to ID periodicity in 3D for the flow between 
parallel planes and to 2D periodicity in 3D for the flow in a rectangular channel) taken at 
two different spatial points [3U] • Another important class of problems associated with the 
quasi-periodic solutions of Laplace equation arises in various problems in electrostatics 
and elastostatics EH EE] ■ 

The relevant representations of Goa for the Laplace equation can in principle be ob- 
tained by taking the limit a — ► in the resulting expressions for the Helmholtz equation. 
In 3D, Hq (cr\r — r' + r s |) in the Schlomilch series (£[]) exhibits a regular limit 

h { n\a\r-r' + r s \) -> —. \ r asa^O (121) 

u v 1 17 |r-r' + r s | v ' 

[see Eq. (j23)] and corresponding Gq (a, R) goes smoothly to the free-space Green's func- 
tion of 3D Laplace equation. However, Hq (o"|r— r'+r s |) displays a logarithmic singularity 
in the same limit (see Eqs. (9.1.3), (9.1.13) of Ref. [HI)- Hence G^(a, R) does not reduce 
to the free-space Green's function Gq(R) = (l/27r)lni? of 2D Laplace equation. Sur- 
prisingly enough, in the case of dual and Ewald representations of Goa the logarithmic 
singularities cooperate in such a way that the limit a — > turns out to be regular even 
for d = 2 (see below). 

In the following, the limit o — > will be discussed in the case of dual and Ewald 
representations of Goa, and in the case of lattice sums Doo in 3D. In taking the limit, 
both Goa and will be considered formally as functions of two independent variables 
a and ky. A reason for doing so is, for instance, solving an implicit equation f|12j) . The 
limit k|| — > will be, if possible, considered afterwords. 

Note that the case of ID periodicity in the x-direction in 2D is rather academic in 
what follows, since in the latter case the Green's function can be calculated in a closed 
form [^U] 

I 

2 



G 0A (R) = -ln 

7T 



7T 



sin -(x -My) . (122) 
a 
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6.1 Dual representations 

Since K± n — > z|kii + k n | in the limit a — > [cf. Eq. 1)170] . in all quasi-periodic cases 
the limit is established [see Eq. flgBJ] by replacing 7i£(K± n \R±\) in Eqs. (j2H), (ESI) 
with Wq (i|k|i + k„||i?j_|). However, for purely imaginary argument ix with x > the 
Hankel functions 7Y^(ix) are related to modified Bessel functions of third kind (see Eqs. 
(9.6.4) and (10.2.15) of Ref. JH|). This results in rapidly decaying terms and exponential 
convergence. As it has been alluded to earlier, absolute exponential convergence of the 
dual representations can only be established under the assumption of R± ^ 0. 

However, the resulting dual representations are singular in the limit kn — > 0. Then 
K±n ~^ ik n and the denominator in Eqs. ()2~Tj). P3*|) vanishes for k„ = 0. 

6.2 Ewald representations 

The a — > limit can also be easily taken in the respective Ewald integral representations 
[)38]). (HOJ), ()4*T|) of free-space quasi-periodic Greens functions: simply substitute in the 
above expressions <p n = and K\ n = — |ki + k n | 2 . Using that the integrals in the series 
over A can be expressed in the a — > limit in terms of incomplete gamma functions (Eq. 
(6.5.3) of Ref. [16J), one finds for the respective Ewald representations (JHHJ), ()30|) . flS]) 
the following expressions: 

G 0A (k|,,R) = — V e 4 ( k n +k "> R ii / C -i e ^(l k n+ k "l 2 C+l^l 2 /C) rfC 

1 /7T\V2 



4tt 2 V2 



r s 6A 



v / 2r[l/2, 


|R + r. 


,•17(2??)] 


1] 


R + r s | 





^ e -y s v^lV^^ri s | LA_iZj (123) 



for ID in 3D, 



1 //7T\ 1 /2 /*°° 

G 0A (k||,R) = — -) V e^i+ k "> R ii / r V2 e -i(|k ll+ k„|^ + |fl x |VC) dC 

^ c -ik rr ., V2ni/2,lR + r s l 2 /(2^)] (m) 



k„gA* 

1 /7T\l/2 



4tt 2 V2 



r s £A 



for 2D in 3D and 
GW(k|i,R) = 



I , x i /o /*oo 

( ) e i(k ii +k ' l) - R ii / ^-i/2 e -i(|k||+k„| 2 c+|i?xl 2 /C)^ 

2 k„eA* J» 



2-kvo \2. 

l^e- k ii-r[0,|R + r s | 2 /(2r / )] (125) 



47T 

r s eA 
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for ID in 2D. The Laplace limit of the respective Ewald integral representations then 
follows straightforwardly by letting ky — > in the above expressions ()123|) . (|124|) . ()125|) . 

Regarding convergence speed, in virtue of the asymptotic T(a, z) ~ z a ~ 1 e~ z for z — > oo, 
|arg z\ < 3tt/2 (Eq. (6.5.32) of Ref. [T^]) the respective Ewald representations remain to 
be exponentially convergent in the a — > 0, ky — ► limit. Note in passing that T[l/2, |R + 
r s | 2 /(2?7)] in Eqs. (|123|) . (|124|) can be expressed via error function [Eq. ()115|) ] as 

r[l/2, |R + r s | 2 /(277)] = v^erfc(|R + r s |/^) (126) 

(Eq. (6.5.17) of Ref. PH]). An alternative exponentially convergent series for G a for 
ID periodicity in 3D in the Laplace case has also been obtained earlier by Linton (see 
series in Eq. (3.26) of Ref. 30 ). However, our expressions have been derived without any 
artificial regularization in the form of a convergence ensuring logarithmically divergent 
series (cf. Ref. [H]). 

Additionally, absolute exponential convergence of the respective Ewald representations 
can also be established for R± = 0, which for numerous alternative representations of G*oa 
in fluid dynamics provides a problem [201 ■ Again, the case R± = can only be attained 
when R)| ^ A. Otherwise the respective Ewald representations become singular. In 3D 
this singularity is explicit, since for some r s G A the denominator |R + r s | vanishes. For 
ID in 2D one has T[0, |R + r s | 2 /(2?7)] — > T(0) for some r s £ A, and the singular behavior 
follows from the pole of the gamma function T(z) for z = 0, or more precisely from the 
asymptotic (Eqs. (6.5.15), (5.1.11) of Ref. 

00 (—z) n 

r(0,z) = E 1 (z) ~ -7-lnz- '— asz-^0, (127) 

n=l 

where 7 is the Euler constant (Eq. (6.1.3) of Ref. [TS] ) . 

As a final remark, note in passing that each of the above Ewald integral representa- 
tions ()123|h (|124|) . ()125|) can be regarded as a one-parametric continuous spectrum of the 
representations for G a- The corresponding dual representations of Sec. 16.11 for kn 7^ 
can be then recovered in the limit 77 — > 0. Indeed, using Hobson's integral representation 
()139|) . the integrals in the series over A* can be expressed in the limit 77 — > in terms 
of the modified Bessel functions of the third kind /Co {K for d — = 2 and Ki/ 2 for 
d — c?a = 1 ) . 
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6.3 Lattice sums in 3D for / = 



In the case of lattice sums Dl, they are defined as expansion coefficients of free-space 
quasi-periodic Greens functions in terms of regular cylindrical (in 2D or spherical (in 3D) 
waves Ji(aR) 3^l(R) [see Eq. (JHJ) above)]. Since unless I = one has Ji(aR) — > in 
the limit a — > [see Eq. ([109)) ]. the lattice sums Dl become singular in the limit for 
1^0. In the case of I = it is reminded here that Dqq (7 function of Karpeshina [T3*l l47p 
determines the energy operator in the one-particle theory of periodic point (zero range) 
interactions and that the Dqq determines the underlying spectrum according to Eq. (|12|) . 

From a mathematical point of view, upon substituting (}2"4*j) for Hq in (fTTj) and assum- 
ing for a while an integer lattice, the lattice sum Dl for / = in 3D can be expressed in 
the limit a —>■ via Epstein zeta function E] 



I, 



3 iki|-r„ 



r„£A 



as 



A) (0,k||) = 
and [Eq. 




k 



(x,i), 

-> as a 



(128) 
(129) 



where we have used Eq. (jT44j) and [Eq. (fT45j) ] that C — ► as a — >■ in 3D. In Eq. 

for g?a = 2, x is a positive quadratic form defined by the scalar product of the basis vector 

of A, Xij — e i " e j- (F° r d\ — 1, all the quadratic forms x reduce to an absolute value). 

The Epstein zeta function in Eq. ()128j) converges absolutely for Re v > d\ and can be 

analytically continued to an entire function in the complex variable v unless ky e A*, in 

which case the Epstein zeta function possesses a simple pole for v = d\. 

Interestingly enough, for dx — 2, d — 3, v — 1 the Epstein zeta function can be 

expressed in a closed form in terms of Jacobi theta functions jT^j. For a ID periodicity in 

3D with a period a one then, starting from Eq. (jH|) , obtains in the limit a — > 

1 



Ax)(0,k,| 



47ra 



In [2 - 2cos(k||a)] . 



(130) 



The expression is obviously singular for kii G A*, in accordance with the singularity of the 
Epstein zeta function for v — d^ — 1. 

Regarding the Epstein zeta functions, note that one could have written dual represen- 
tations in the limit a —>■ for R± = 0, d = 3, and provided that kn A*, as 



GoA^knR) 



e ik r R|j 
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R 



(x,i), 



(131) 



where 



Z 



k|| 



k„eA* 



k|| + k n 



v 



(132) 



7 Discussion 



7.1 The choice of the Ewald parameter r) 



Each of the exponentially convergent Ewald representations (J3*8|) . (|4"0|). (jH)) for G a can 
be viewed as a one-parametric family of representations. A corresponding image-like se- 
ries (JHJ) and a dual representation [Eqs. (|2"Tj) . (|2*HJ) ] can be seen then as two ends of the 
one-parametric continuous spectrum of the representations for Goa- Obviously, by varying 
the point rj, at which the integration is split, the convergence characteristics of the repre- 
sentation can be altered. In most cases, the value of the Ewald parameter 77 is chosen to 
balance the convergence of respective and contributions. This leads occasionally 
to the criticism that an arbitrary optimization parameter enters the evaluation of lattice 
sums. On contrary, the invariance of D^s on the value of Ewald parameter 77 serves as a 
check of a correct numerical implementation. The Ewald parameter 77 can often be varied 
by several orders of magnitude without affecting the results in a wide frequency window. 
However, for some range of 77 values one can enter a numerically unstable region: the 
respective and contributions have opposite signs and similar magnitude, which 
is several orders larger than the magnitude of resultant D^. This instability can easily be 
remedied by the choice of some other value of rj, or one can follow the recipe of Berry [21] 
and chose 77 to depend on a and I, and thereby prevent numerical instability completely. 
Indeed, although the results presented here have been obtained by a uniform /-independent 
choice of 77, one can easily modify the above derivations to the case of /-dependent r\ |24j . 

7.2 Numerical convergence 

Exponentially convergent representations of lattice sums summarized here provide a sig- 
nificant advantage in terms of computational speed, while maintaining accuracy, over 
alternative expressions of lattice sums. In the special case of a ID periodicity in 2D this 
is demonstrated in Table I below. In the case of ID periodicity in 2D, even with the 
latest progress due to Yasumoto and Yoshitomi [33], it took 40 seconds to compute G a 
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on SPARC workstation from lattice sums with 14 digits accuracy at a single point and 
frequency. This was in striking contrast to the calculation of exponentially convergent 
lattice sums in the so-called bulk cases, i.e., when Goa is periodic in all space dimensions. 
The lattice sums for an infinite 2D lattice in 2D |23j and an infinite 3D lattice in 3D [2|. 
The respective convergence times (on a PC with Pentium II processor) for a set of bulk 
2D and 3D lattice sums with 6 digits accuracy are less than « 0.03 second (for a cut-off 
value of l max = 20) and ~ 0.8 second (for a cut-off value of l max = 6) [BT] . 

TABLE I. Quasi-periodic Green's function Goa for off-axis incidence at an angle 6 = tt/8 
upon a ID lattice oriented along the x-axis in 2D with X/v = 0.23. Here Vq is the length 
of a period (primitive lattice cell) along the x-axis. In rows labeled by D, values of Goa 
obtained by a direct summation of its dual representation (spectral domain form) (taken 
from Tab. Ill of Ref. [S3])- These data are compared against those in rows labeled by E, 
obtained by the (complete) Ewald-Kambe summation [Eqs. (jHEJ), ()103j) . ()119j) with the 
Ewald parameter rj = 0.011]. Data in the respective rows labeled by YY and NMcP are 
those obtained by Yasumoto and Yoshitomi and Nicorovici and McPhedran |3*3*j 

methods. 





X 


y 


Re G a 


Im G 0A 


D 
E 

YY 
NMcP 


0.2 


0.03 


0.117120006144932 
0.117120006144932 
0.117120006144941 
0.117120006141860 


-0.108131857633201 
-0.108131857633206 
-0.108131857633206 
-0.108131857633197 


D 
E 

YY 
NMcP 


0.2 


0.003 


0.115891895634567 
0.115891895634565 
0.115891895634577 
0.115891895630095 


-0.103497063599642 
-0.103497063599651 
-0.103497063599646 
-0.103497063599643 


D 
E 

YY 
NMcP 


0.2 


0.0003 


0.115881138140449 
0.115881138140448 
0.115881138140457 
0.115881138135960 


-0.103450147416784 
-0.103450147416794 
-0.103450147416788 
-0.103450147416785 



The computational time to reproduce a value of Goa in Table I with accuracy of within 
8 x 10~ 15 of that obtained by a direct summation turns out to be ~ 0.2s, in line with 
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the respective ~ 0.03s and ~ 0.8s for convergence time of a set of bulk 2D [o2| and 3D 
lattice sums [HI] with 6 digits accuracy. This should be compared to 1232s of Nicorovici 
and McPhedran [33J, or, to 40s of Yasumoto and Yoshitomi [33] (the computational times 
have been taken from Ref. 35J). The exponentially convergent representation [Eqs. ()86j) . 
f!103|) . ([119)1 ] (i) can be implemented numerically more simply and (ii) converges roughly 
200 times faster than the previous best representation 35J. Of the cases tested, the 
simplest case of a constant rj = 0.011 was chosen. 

For a 2D periodicity in 3D, a comparison of the speed and accuracy of exponentially 
convergent representation of lattice sums [Eqs. (J85|) . (jl02|) . ()118|) ] with respect to al- 
ternative expressions for Goa has been summarized in Ref. [2Zj. Again, exponentially 
convergent representation of lattice sums turns out to be convergent for a given accuracy 
significantly faster. 

The reader is invited to perform some additional tests by using several publicly avail- 
able F77 codes. In the case of a 2D periodicity in 3D, numerical codes can be obtained 
from Comp. Phys. Comp.: for a complex 2D lattice in 3D see routines DLMNEW 
and DLMSET of Ref. [HZ]; for a simple 2D Bravais lattice in 3D see routine XMAT 
of Ref. |72] • The above codes have been implemented in electronic, acoustic, and 
electromagnetic LKKR codes and successfully tested time and again in various cases 
|221in3inniin71inHll7nil721l7S]- A limited Windows executable which incorporates the 
lattice sums within a photonic LKKR code for the calculation of reflection, transmission, 
and absorption of an electromagnetic plane wave incident on a square array of finite length 
cylinders arranged on a homogeneous slab of finite thickness is available following the link 
http:// www.wave-scattering.com/caxsrefl.exe 



F77 source code for a ID Bravais periodicity in 2D is freely available at 
http:// www.wave-scattering.com/ola.fl (implementation instruction are described on 



http:// www.wave-scattering.com/dlsumlin2.html). The code has been implemented in 



a corresponding LKKR code and successfully tested against experiment in Ref. |7ij . 
A limited Windows executable calculating the reflection, transmission, and absorption of 
an electromagnetic plane wave incident on a square array of infinite length cylinders in the 



plane normal to the cylinder axis, is available following the link http: / /www. wave-scattering.com/rtalin21 
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7.3 Outlook 



The present work can be straightforwardly extended in several directions. First, as in bulk 
cases [H |2], for a 2D periodicity in 3D and for a ID periodicity in 2D [3], the condition 
of a simple lattice for a ID periodicity in 3D can easily be relaxed to an arbitrary periodic 
lattice. Note that the case of a non-Bravais lattice additionally requires the calculation of 
the series (0) with the origin of coordinates displaced from the lattice by a fixed non-zero 
vector. Consequently, the term involving r n = is no longer singular. Therefore, in 
the latter case the lattice sums are expressed as the sum of solely and D^, where 

(2) 

D K L does include the r n = term. These supplementary series can easily be determined 
following the recipes of Refs. [H El El H] • 

Second, following the work of Ohtaka [68J and Modinos [69 j for a 2D periodicity in 
3D, in the vector case of electromagnetic waves for a ID periodicity in 3D (2D case 
is trivial as it reduces to a scalar problem), the lattice sums and structure constants 
can easily be obtained from those in the scalar case presented here. It is only required 
to multiply the scalar quantities with appropriate numerical factors of geometric origin 
[SHI EDI IHH1 EH]- This possibility is a consequence of a fact, as first shown by Stein JED], 
that vector translational addition coefficients can be derived from pertinent scalar addition 
coefficients. A more involved, but possible, is a generalization of the presented results to 
semi-infinite cases, when periodicity is imposed on a half-line, or in a half-space only 

El EH EH- 

8 Summary and conclusions 

A classical problem of free-space Green's functions Goa representations of the Helmholtz 
equation was studied in various quasi-periodic cases, i.e., when an underlying periodicity 
is imposed in less dimensions than is the dimension of an embedding space. Exponentially 
convergent series for the free-space quasi-periodic Goa and for the expansion coefficients 
Dl of Goa in the basis of regular (cylindrical in two dimensions and spherical in three 
dimension (3D)) waves, or lattice sums, were reviewed and new results for the case of a 
one-dimensional (ID) periodicity in 3D were derived. The derivation of relevant results 
highlighted the common part which is applicable to any of quasi-periodic cases. 
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Exponentially convergent Ewald representations ([SEj) , (001), PTj) for G a (see also Sec. 
16. 2|) and for lattice sums Dl hold for any value of the Bloch momentum and allow Gq\ 
to be efficiently evaluated also in the periodicity plane. After substituting the resulting 
expressions for [Eqs. 1)55). ljB5). ©], ^S? [Eqs. dEJ, and L>[ 3) [Eqs. (fTTHj) . 

()118p . ()119|) ] into defining equation (jlS)) for .D^, an alternative exponentially convergent 
representations for Schlomilch series (|TT|) of cylindrical and spherical Hankel functions of 
any integer order are obtained. 

The quasi-periodic Green's functions of the Laplace equation were studied as the lim- 
iting case of the corresponding Ewald representations of G \ of the Helmholtz equation 
by taking the limit of the wave vector magnitude going to zero. Thereby, exponentially 
convergent representations of Goa i n the Laplace case were obtained, which are conver- 
gent (unless R G A) also in the periodicity plane. An alternative exponentially convergent 
series for Goa for ID periodicity in 3D in the Laplace case has also been obtained earlier 
by Linton (see series in Eq. (3.26) of Ref. [30 ). However, our expressions have been 
derived without any artificial regularization using a convergence ensuring logarithmically 
divergent series (cf. Ref. .30J). 

The results obtained can be useful for numerical solution of boundary integral equa- 
tions for potential flows in fluid mechanics, remote sensing of periodic surfaces, periodic 
gratings, in many contexts of simulating systems of charged particles, in molecular dy- 
namics, for solving the spectrum of particular open resonators, for the description of 
quasi-periodic arrays of point interactions in quantum mechanics, linear chains of spheres 
and nanoparticles in optics and electromagnetics, and of infinite arrays of resonators cou- 
pled to a waveguide, and in various ab-initio first-principle multiple-scattering theories for 
the analysis of diffraction of classical and quantum waves. 
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A Integral representations of 

For the Ewald summation with r s = excluded, one uses the Schlafli integral representa- 
tion, 

Hl 1 \z) = — [ u-^e^^-^du, (133) 
m Jc- 

where the contour C_ is the contour which goes from the origin to 5 > 0, continues along 
a semicircle in the upper half-plane with radius 5 to —5, and goes along the negative 
real axis to infinity (see Fig. This results in an exponentially decreasing integrand at 
the integration contour ends for Re z > 0. The Schlafli representation, which yields the 
Hankel functions Hu (z) as moments of the generating function of the Bessel functions 
Ji (z) (see [7j, p. 14), is obtained upon substitution u = e* in the integral representation 
(see (9.1.25) of Ref. [E|), 

-i /*oo+7n 

H^(z) = — e zshlht ~ ut dt, (134) 

m J-oo 

taken along the contour shown in Fig. HJ which is valid for | arg^| < tt/2. 
Setting in the Schlafli representation ()133j) z = or and 

2t / /t \ -l 

" e, du=-2[—) ^, (135) 



a * ' V2r 
one arrives at 

HP {or) = 2 ±^(1Y f r 2,-i e -r^ + .W) dC (136) 
7r V2r/ Jq 

where C is the so-called Ewald contour (see Fig. |3J), which leaves the origin along the ray 
arg £ = arg a — it/A, then returns to the real axis and continues along the positive real 
axis to infinity. 

Upon substituting in the Schlafli representation (|133|) z = ar and 

u = -(, (137) 
a 

one arrives at 

H<»(ar) = -.{-)" I C v - l e^ C -^ /0 dC (138) 

Let us consider for a while a as a general real parameter. Then, following discussion at 
the end of Appendix 2 of Ref. j2I], the Schlafli integration contour for a positive a 2 can 
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be deformed to that shown in Fig. |21 Afterward, with the use of Jordan's lemma 
p. 115) integration contour can be deformed to that from to ioo on the imaginary axis. 
For a negative a 2 one would then, as the result of the substitution ()137|) . arrive at the 
contour from to oo along the positive real axis (see Appendix 2 of Ref . [21] ) . 

Eventually we provide Hobson's representation for the modified Bessel functions of the 
third kind [H] 



K v {kq) = I" C~ l e-^^ WK) dC (139) 



B Properties of harmonics 34 

Throughout this paper complex harmonics (cylindrical, Yj = e d ^/ \/2~k, for d = 2 and 
spherical for d = 3 ^01 IH 1221 EH1 ; for ID harmonics see Ref. j5U]) are used. Under complex 
conjugation, they behave according to 

r yi(k), id 

yi(R)=l y_,(R), 2D (140) 
I (-l) m ^- m (R), 3D. 
In 3D the property goes under the name of the Condon-Shortley convention. In any 
dimension, the harmonics satisfy the inversion formula 



the orthonormality 



and closure 



y L (-r) = (-l) l y L (i), (141) 
j>y L {x)yi,{r)dSl = 8 LL ,, (142) 
^^L(r)K(r')=fe(r,r'). (143) 



Here 5n(r, r') is the delta function on the unit sphere whereas dQ is the usual angular 
measure which is determined by the relation dr = r d ~ l drdQ. (In ID the angular integral 
§ dQ reduces to the summation over the forward and backward directions.) 
Note that in any dimension, 

yo = ^=- (144) 
VA 



Since y is a constant, combining Eq. ()144|) with the orthonormality ()142|) of yields 

j>y L {v)dn = Vas lo . (145) 



41 



B.l The plane- wave expansion 

An expansion of the plane-wave expansion in the angular momentum basis is 



e*- = A Y,^J\l\{kr)y L {i)ym = Aj2^Jdkr)yi(i)y L (k), (146) 

L L 

where 

A= (pdQ= { 2tt, 2D (147) 




where dfl here is the usual angular measure. The series ()146|) converges uniformly as |k| 
and r run through compact sets of 1R and 1R 3 (Theorem XI.64f of Ref. [ZH])- The absolute 
value of I is used in f!146|) in case the sum over angular momenta in 2D runs from minus 
to plus infinity. 

The plane-wave expansion ()146|) is also valid for complex arguments. It is interesting 
to note that yi is no longer the complex conjugate of for complex k (e.g., in 3D 
because of the complex nature of associated Legendre functions). However, the relations 
()140|) remain the same as in the case of harmonics of a real argument jU ESI- For 3D case 
see, for instance, Appendix 1 of Ref. jl] or Appendix A of review by Tong [66J. The 2D 
case follows by a straightforward adaptation of the 3D case, whereas the ID case is trivial 

C Free-space scattering Green's function 

One has [see Eq.fJHJ)] 

G >,r,r') = JL/ J^l ie d^^-iCHt{a\v-v'\l (148) 



where 



7T 



A 2o-' 



^ ID 



C = ^a"=( I 2D (149) 



2 (2ir) d 



is a real positive number for positive energies. The 2nd equality in (jSJ) is established 
by expanding the exponential into regular waves according to Eq. ()146|) . performing the 
angular integral using the identity (|145|) satisfied by the harmonics 3^, and eventually 
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performing the remaining radial integral using the integral identity 

' = [J^^'^-^XH, (150) 

which can be easily established by contour integration in complex plane. In arriving at 
Eq. (jl50j) one substitutes for J according to 

Ji(kr) = l[H+(kr)+H7(kr)]. (151) 

and applies the identity (analyticity property) 

H^ze'™) = (-l) l+d+1 Hf(z) (152) 

(see Eqs. (9.1.16) and (10.1.18) of Ref. [S] and Ref. [SU]) for I = 0. 
One can show that 

|x-y| ~ Ixl - ^^ + C(l/|x| 2 ) aslxl^oo. (153) 

|x| ' 

Therefore, upon using the asymptotic properties (0) of Hq (<t|x — y|) for |x| — > oo, 

G+(x,y) ~ ^(IxDe-^/W = / CT (|x|) e - ik '^ (154) 

where 

/ CT (|x|) ~ G + (x,0)| x ^ oo = -iCH+(v\x\), k' = W|x|. (155) 

Function / CT (|x|) describes outgoing waves in a given dimension. Explicitly 

-^e iCT l x l, ID 
/.(|x|)= "v^— 2D (156) 



1 elCT|x| 3D. 



4tt |x| ' 

The product AC in the partial wave expansion (JHJ) of the free Green's function can be 
independently determined from the condition that the discontinuity of radial derivatives 
of the free Green's function at coinciding arguments multiplied by r d ~ l is exactly one. 
The factor r d ~ l follows from a fact that the integral measure dr can be written as dr = 
r d ~ 1 drdQ. Since the harmonics are orthonormal in the measure dfl [Eq. ()142j) ]. 
the discontinuity of radial derivatives of the free Green's function can be conveniently 
expressed by the Wronskian 

W x [f(ax),g(ax)] = f(ax)g'(ax) - f(ax)g(ax), (157) 
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where prime denotes first derivative with respect to x. Now 
for z = ar and 

W r [J l {ar),nl{ar)} = iW r [Ji(ar),M(<^))- 
Knowing the Wronskian 

f 1, ID 

W g [J t (z),M(z)] = JiNl - JIM = I 2/(nz), 2D (158) 

{ 3£>, 

the product AC can be then found directly from the discontinuity given by the relation 

AC = {ar d - x W z \Ji{z) Ml{z)\}~ X . (159) 

D Jacobi identities 

The equivalence of the lattice sum [see Eq. (jHJ)] and the eigenvalue expansion [see Eq. 
1)13)1 ] in the case of the heat equation with the Bloch boundary conditions in a box leads 
to the identity 

K A (R,t) = — i— - V e *-r. e -(R-r.) a /(«) = - y e -(k+k n )V(k + k„).x_ (160) 

In a special case, for a simple cubic lattice with a unit lattice constant one has A = A*. 
Upon substituting k„ = 2im, n being an integer valued vector, and R = k = 0, identity 
(fTO yields 

1 Ee -^., = Ee -w« (161) 

^ ' neA neA 

For d = 1, the latter is a special 6 = case of the famous number theoretical Jacobi theta 
function identity (upon rescaling t — > 47r£), 

oo oo 

exp(-7rn 2 t - 27rm0) = r 1/2 ^ exp [-tt(Z + #) 2 /t] , (162) 

n=— oo !=— oo 

which is valid for a complex 9 and Re £ > 0. The Jacobi theta function identity can 
also be proved by applying the Poisson sum rule and is also sometimes referred to as 
Poisson- Jacobi formula. 
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An alternative form (|3"5Jl of the Jacobi formula, as has been used by Kambe |4*1 |2T)1 12~T] . 
is obtained upon substitution t = (/2 as is valid for Re ( > 0. After the substitution 
t = l/(4£ 2 ), the generalized Jacobi identity (jl60j) yields 

^2 e -(R-r s ) 2 £ 2 +ik.r s _ ^_ ^ e ~ (k ' l+k)2/(4?2)+i(k+k " ) - R . (163) 

r s 6A k n £A* 

In like manner, upon the substitution R — > — R, k — > — k, k n — > — k n , one finds 

g-(R+r s ) 2 e 2 -ik.(x+r s ) _ ^_ e -(k„+k) 2 /(4? 2 )+ 4 k rl -R > 
r s eA k„eA* 

The last two identities are often referred to as the Ewald identities (cf. Refs. CHI QSj)- 

E McRae derivation of dual representations 

In the quasi-periodic case, the dual representation can also be established by applying 
an Ewald integral representation of Green's function and the generalized Jacobi identity 
This path has been originally followed by McRae [49J for a 2D periodicity in 3D. Here it 
will also be outlined for the remaining two cases. 
Let 

S djdA (a,r) = ^H+(a|r + r s |) e - k -^ +r ), (165) 

r s eA 

where A is as usual a G?A-dimensional lattice (d\ < d) and r s are the lattice points. Then 
[cf. Eq. ©] 

G 0A (a,r,r') = - i ^^- d a d - 2 S d , dA (a,r-r'). (166) 
Let us first consider a special case 

Ss, dA = £ h§-\a\r + r s \)e-^ +r l (167) 

r s eA 

Now the Ewald integral representation [Eq. (J136)) of Appendix IX] 

h^(ar) = -^L [ e ^/(^)-^ d£ (168) 
V 1Tcr Jc 

is substituted for in Eq. ([167)1 . Here the path of integration leaves the origin along 
the ray arg £ = arg a — tt/4 and then returns to the real axis. Afterward one finds 

S, dA = -JL y [ e ^/(« a )-[(r.+r„)' + ri 1^-^+0 ^ (1 g g) 
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By regarding the integral as the limit Im a —>■ 0+, the double summation here can be 
shown to be absolutely convergent for r± ^ 0, and hence the order of contour integration 
and summation can be reversed. Indeed, on the interval (rj, oo), rj > 0, along the positive 
real axis the sum of integrands is less than 

const x e" r ^ 2 • (170) 

On the other hand, in the proximity of zero along the Ewald integration contour £ = 
e~ m 'H, where t £ (0, rj) is a real positive number, and hence the sum of integrands is less 
than 

x e- CT "/* 2 , (171) 



const _ T » .■/■-: 

V P 

1 it 2 

1 — e u 



where a" denotes the imaginary part of a. After the reversal of summation and integration, 
the respective Jacobi identities ()164|) for d\ — 1 and d\ = 2 can be applied. To this end, 
the respective expressions for d\ = 1 and d\ = 2 are formally identical and they only 
differ in the dimensionality of the lattice A. 

Upon using the generalized Jacobi identity ()164J) for d\ = 1 with R = v\\ one obtains 

S 31 = -— y e^ r n / e K U^)-r\e (172) 
^kfil* Jo 

where K± n is as usual, given by Eq. fjlTj) . and A* is a corresponding reciprocal ID lattice 
with lattice points k n . Now, for r > and | argcr| < n/2, Eq. ()136j) yields 

if £V) = - [ e^ 2 ^ 2 T 1 ^. (173) 



c 



Therefore, 

7T 

k n eA* 

For g?a = 2, upon using the generalized Jacobi identity (|164j) with x = rii, one obtains 



S 3 .i = — V e ik »- r Hi7 (1) (^ ±n |r ± |). (174) 

WOO" ^ 



fta= Z!![ y e *»* 2 /" y e ^i n /(4^)-rie 2 r 2^_ (175) 



Upon substituting 

r 2 = / e-^dt, (176) 



00 

2+ 



( I 
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which is a special case of the Euler transformation ^UJ , interchanging the order of 
integration and using [Eq. ()168|) : Appendix IX] 



2 

7T 



e «*m*)-r>e dfj (177) 

C 



the integral on the r.h.s. of Eq. (jl75j) becomes 

r r°° ciK±n(r 2 , +i) 1/2 



/7T 

Upon substituting 



j/' M '^L (l78) 



u = (ri + 1) 1/2 , du = dt/[2{r 2 x + 1) 1/2 }, (179) 



one eventually arrives at 

roc iK^irl+t) 1 / 2 roo JK ±n \r±\ 

/ 7 2 — m~2 dt = 2 / e iK±nU du = 2i— = -2\r ± \h^(K ±n \r ± \). (180) 

JO V± + J J Jr ± J^±n 

For real K± n one adopts as the value of the integral the limit as Im K± n — ► 0+. Eventually, 

^3,2 = — V e <k --'ii|r JL |/i5 1) (JifjL„|rj.|). (181) 

For completeness, and without a proof, we also present the result for 

S 2A = V H^\a\r + r s |)e- lk ^ +r ) = - V e^n |r J^i^r J), (182) 

— fn — 

r s eA k n eA* 

which have been used as an intermediary step in the derivation of results of Ref . j2E] ■ 

When the respective dual representations (|174|) . (|181|) . (|182|) are substituted into Eq. 
(|166j) for a corresponding G a, r in Eqs. (|174j) . (|181|) . (|182|) is substituted by R. 

F General properties of free-space quasi-periodic Greens 
functions and of lattice sums 

For any r s G A, k n G A*, Gq\ satisfies the following trivial properties, 

G 0A (a,k,|,R) = G 0A ( C r,k||,R + r s ) = G 0A (a,k|| +k„,R). (183) 

Obviously, Goa is only a function of the projection kii of k upon A* (note that kii = k for 
c?a = d). 
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Except for the frequencies which satisfy a 2 = (ky +k n ) 2 for some k n G A*, Goa satisfies 
the following reflection symmetry property [2H E] , 

G 0A ((T,k||,R||+R ± ) = Gf 0A ( < 7,k||,R|| -R ± ). (184) 

From dual representations (J3TJ), (f^Hj) it follows that the respective Goa are hermitian 
for the interchange of variables ry and rji and complex symmetric for the interchange of 
variables r_|_ and r' ± . 

Obviously, for any r n G A also — r n G A. Upon combing the inversion formula (jl41|) 
of angular- momentum harmonics 3^ with the defining Eq. (fTTj) for Dl one finds 

r„eA 

The summation here is performed over the subset A C A of equivalence classes of r n G A 
with respect to spatial inversion r n — > —r n . Thus, for I even the corresponding lattice 
sums Dl(o~, ky) are even functions of the Bloch vector ky, 

D L {a,k {l ) = D L {a,-k l{ ). (186) 

In the special case of a 2D periodicity in 3D one has additionally 

D L (a,k {l ) = D L (a,kf) = D L (a,kJ), (187) 

where the respective wave vectors kjj" = (k 1 ,—k 2 ), kjj 1 = (—k 1 ,k 2 ) are formed from the 
Bloch vector ky components k± and k 2 . 

In the special case of Doo(a, ky) for a = iy/—z there are known some further general 
analytic properties in the complex z-plane, which have been established by Karpeshina 
[13 EI: 

• Doo(i^— z, ky) is for a fixed Bloch vector ky an analytic function on the complex 
plane cut along the real half- axis z > |ky| 2 with Im D Q0 ^ on both sides of the 
cut. 

• On the real half-axis z < |k||| 2 the function D^i^J— z, ky) is real, smooth, and 
increases monotonically (0 < d z D 00 (i\/—z, ky) < oo) from — oo to oo. 
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G Alternative definitions of lattice sums and struc- 
ture constants 

If instead of D\ of Eq. (J2J) the difference 

V A (a,k h R) =GoA(a,k,R)-G (cx,R), (188) 
which is also regular for R — > 0, is expanded in terms of the regular waves, 

D A (<7,k,|,R) = Y, V L(°^\W aR )yd&) (189) 

L 

= Ys9LL'^M)^ r )M°r')y L (i)yUi'). (190) 

L,L' 

this results in alternative lattice sums T>l and structure constants Qlv ■ We recall here 
that the structure constants can in a known way be unambiguously determined from the 
lattice sums |5lj . 

The lattice sums T>l and Dl and structure constants Atv and giu, where 

D A (a,k h R) = Y,A L L\°M)^ r )M°r')y L {i)yW)i (191) 

L,L> 

are related to each other as follows: 

V L (a, k„) = D L (a, k„) + iCA l ' 2 5 L0} (192) 

g LU (a, k,| ) = A LU (a, k|, ) + iAC 5 LU . (193) 

A and C here are the familiar numerical constants which have been defined by Eqs. f|147|) 
and ()149|) . respectively. Relations ()192J) and ()193|) follow easily from 

iImG+((7,R) = -iCJ (aR) = -iCA 1 ' 2 J (aR)y {R) 

= Y.9LLiJi(°r)Mori)y L {T)yi,{r'), (194) 

LL' 

where 

g° LL , = -iAC5 LLI . (195) 

In going from the first to the second equality in (jl94j) we have used that in any space 
dimension y (R) = A~ 1/2 [Eq. (ITU) ]. The final expressions then readily follows from the 
partial wave expansion of the free Green's function (see Eq. Q of Appendix IUj) . 
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Figure captions 



Figure 1 - Geometry and parameters - A plane wave with wave vector k incident on A 
with an incidence angle 9. 

Figure 2 -a) Deformation of the integration contour for K\ n > in Eqs. (J3*8|) . (j4TJj) . (}4*T]) 

before Jordan's lemma is applied to the quarter-circles, b) For Kj_ n < the integration 
contour in Eqs. (jSSj), pUj). (}4*T]) can be considered as a sum of two contours, one from r\ 
to zero and the other from zero to +oo. 

Figure 3 - Integration contour in the Schlafli integral representation of 

Figure 4 - Rectilinear integration contour used in the derivation of the Schlafli integral 

representation of Hv . 

Figure 5 - Integration contour in the Ewald integral representation of . 
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